This notebook provides the analysis regarding Hypothesis 1, Prediction 1b for the article: Heyne, M., Derrick, D., and Al-Tamimi, J. (under review). “Native language influence on brass instrument performance: An application of generalized additive mixed models (GAMMs) to midsagittal ultrasound images of the tongue”. Frontiers Research Topic: Models and Theories of Speech Production. Ed. Adamantios Gafos & Pascal van Lieshout.
# specify directory to save models and summaries
output_dir = "updated_models"
# specify whether to run models -> if set to false script will attempt to load saved models from output_dir
run_models = FALSE
load_packages = c("parallel","mgcv", "itsadug", "rlist", "ggplot2", "plotly", "dplyr")
for(pkg in load_packages){
eval(bquote(library(.(pkg))))
if (paste0("package:", pkg) %in% search()){
cat(paste0("Successfully loaded the ", pkg, " package.\n"))
}else{
install.packages(pkg)
eval(bquote(library(.(pkg))))
if (paste0("package:", pkg) %in% search()){
cat(paste0("Successfully loaded the ", pkg, " package.\n"))
}
}
}
Successfully loaded the parallel package.
Successfully loaded the mgcv package.
Successfully loaded the itsadug package.
Successfully loaded the rlist package.
Successfully loaded the ggplot2 package.
Successfully loaded the plotly package.
Successfully loaded the dplyr package.
rm(load_packages, pkg)
# detect number of cores available for model calculations
ncores = detectCores()
cat(paste0("Number of cores available for model calculations set to ", ncores, "."))
Number of cores available for model calculations set to 8.
# plot multiple GAM model outputs
plotly_model_outputs <- function(model, changing_cond, changing_var, constant_cond1, constant_var1, values, constant_cond2=NULL, constant_var2=NULL, print=TRUE){
if(length(constant_var1)>1 | length(constant_var2)>1){
print("Error: Constant variables can only have length 1.")
}else{
if(!is.null(constant_cond2) && !is.null(constant_var2) && length(changing_var)==2){
# works for models Notes.gam...
cond_p1 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[1], "', ", constant_cond1, "='", constant_var1, "', ", constant_cond2, "='", constant_var2, "')")))
p1=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p1)))
cond_p2 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[2], "', ", constant_cond1, "='", constant_var1, "', ", constant_cond2, "='", constant_var2, "')")))
p2=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p2)))
# set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
max_ul = max(p1$fv$ul, p2$fv$ul)
max_fit = max(p1$fv$fit, p2$fv$fit)
maximum=max_fit+((max_ul-max_fit)/2)
# plot in polar coordinates
p=plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=changing_var[1]) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=changing_var[2]) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)),
title=paste0("GAM smooths @", constant_var1, " & ", constant_var2))
p
}else if(is.null(constant_cond2) && is.null(constant_var2) && length(changing_var)==2){
# no specific case yet
cond_p1 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[1], "', ", constant_cond1, "='", constant_var1, "')")))
p1=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p1)))
cond_p2 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[2], "', ", constant_cond1, "='", constant_var1, "')")))
p2=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p2)))
# set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
max_ul = max(p1$fv$ul, p2$fv$ul)
max_fit = max(p1$fv$fit, p2$fv$fit)
maximum=max_fit+((max_ul-max_fit)/2)
# plot in polar coordinates
p=plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=changing_var[1]) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=changing_var[2]) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)),
title=paste0("GAM smooths @", constant_var1, " & ", constant_var2))
p
}else if(is.null(constant_cond2) && is.null(constant_var2) && length(changing_var)>2){
# works for models NZE.gam... or Tongan.gam...
for (i in 1:length(changing_var)){
if(changing_cond == constant_cond1){
cond_p1 = capture.output(cat(paste0("list(", constant_cond1, "='", constant_var1, "')")))
p1=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p1)))
# exception for KIT (='\\\\') when not using IPA symbols
if (changing_var[i]!="\\\\"){
cond_p2 = capture.output(cat(paste0("list(", changing_cond, "='", changing_var[i], "')")))
p2=plot_smooth(x=get(model), view=values, cond = eval(parse(text=cond_p2)))
}else{
p2=plot_smooth(x=get(model), view=values, cond = list(token.ord='\\\\'))
}
# set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
max_ul = max(p1$fv$ul, p2$fv$ul)
max_fit = max(p1$fv$fit, p2$fv$fit)
maximum=max_fit+((max_ul-max_fit)/2)
# plot in polar coordinates
p=plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=constant_var1) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=changing_var[i]) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)),
title=paste0("GAM smooths ",constant_var1," vs ", changing_var[i]))
Sys.sleep(0)
print(p)
}
}
}
}
}
# plot smooths with shading to indicate significant differences (Matthias Heyne, 2019)
plotly_smooths_w_sig_diff <- function(model, condition, var1, var2, values, language, print=TRUE){
# get intervals of significant differences by running plot_diff
# unfortunately setting plot=FALSE doesn't work as intervals of significant difference are not displayed!
# hardcoded condition
# output = capture.output(plot_diff(get(model), view=values,
# comp=list(tokenPooled.ord=c(var1, var2))))
# output = capture.output(plot_diff(get(model), view=values, comp=list(langNoteInt.ord=c(paste0("Tongan.", note, ".", intensity), paste0("NZE.", note, ".", intensity)))))
names_smooths=list()
if (condition=="tokenPooled.ord" && length(language)==1){
output_comp = capture.output(cat(paste0("list(", condition, "=c(var1, var2))")))
names_smooths[1]=var1
names_smooths[2]=var2
plot_title = paste0("GAM smooths ", language, " ", var1, " vs ", var2)
}else if (condition=="langNoteInt.ord" && length(language)==2){
output_comp = capture.output(cat(paste0("list(", condition, "=c('", language[1], ".", var1, ".", var2,
"', '", language[2], ".", var1, ".", var2, "'))")))
names_smooths[1]=paste0(language[1], ".", var1, ".", var2)
names_smooths[2]=paste0(language[2], ".", var1, ".", var2)
plot_title = paste0("GAM smooths ", language[1], ".", var1, ".", var2, " vs ", language[2], ".", var1, ".", var2)
}else if (condition=="native_lg.ord"){
output_comp = capture.output(cat(paste0("list(", condition, "=c(var1, var2))")))
names_smooths[1]=var1
names_smooths[2]=var2
plot_title = paste0("GAM smooths ", language, " ", var1, " vs ", var2)
}
# output_comp = capture.output(cat(paste0("list(", condition, "=c(var1, var2))")))
output = capture.output(plot_diff(get(model), view=values, comp=eval(parse(text=output_comp))))
# no significant difference
if ((length(language)==1 && length(output)==7) | (length(language)==2 && length(output)==6)){
cat(paste0("Smooths for ", var1, " & ", var2, " are not significantly different.\n"))
dat1 = NA
assign(paste0("int_sig_diff_", var1, "_", var2), dat1, envir = .GlobalEnv)
rm(dat1)
# run plot_smooth to grab data for polar plots
p1 = plot_smooth(x=get(model), view=values, cond=list(tokenPooled.ord=var1, tokenPooled.ord=var2))
p2 = plot_smooth(x=get(model), view=values, cond=list(tokenPooled.ord=var2, tokenPooled.ord=var1))
# set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
maximum=max(p1$fv$fit, p2$fv$fit)+((max(p1$fv$ul, p2$fv$ul)-max(p1$fv$fit, p2$fv$fit))/2)
# plot in polar coordinates
p=plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)),
title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
Sys.sleep(0)
print(p)
# there are differences...
}else{
# grab intervals of significant differences from output
if (length(language)==1 && length(output)>=8){
sig_diff1 = c(as.double(unlist(strsplit(unlist(strsplit(output[8], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[8], " "))[3]))
}else if (length(language)==2 && length(output)>=7){
sig_diff1 = c(as.double(unlist(strsplit(unlist(strsplit(output[7], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[7], " "))[3]))
}
if (length(language)==1 && length(output)>=9){
sig_diff2 = c(as.double(unlist(strsplit(unlist(strsplit(output[9], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[9], " "))[3]))
}else if (length(language)==2 && length(output)>=8){
sig_diff2 = c(as.double(unlist(strsplit(unlist(strsplit(output[8], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[8], " "))[3]))
}
if (length(language)==1 && length(output)>=10){
sig_diff3 = c(as.double(unlist(strsplit(unlist(strsplit(output[10], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[10], " "))[3]))
}else if (length(language)==2 && length(output)>=9){
sig_diff3 = c(as.double(unlist(strsplit(unlist(strsplit(output[9], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[9], " "))[3]))
}
if (length(language)==1 && length(output)>=11){
sig_diff4 = c(as.double(unlist(strsplit(unlist(strsplit(output[11], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[11], " "))[3]))
}else if (length(language)==2 && length(output)>=10){
sig_diff4 = c(as.double(unlist(strsplit(unlist(strsplit(output[10], " "))[1], "\t"))[2]), as.double(unlist(strsplit(output[10], " "))[3]))
}
# write intervals of significant difference to variable
if ((length(language)==1 && length(output)>=11) | (length(language)==2 && length(output)>=10)){
dat1 = c(sig_diff1, sig_diff2, sig_diff3, sig_diff4)
}else if ((length(language)==1 && length(output)>=10) | (length(language)==2 && length(output)>=9)){
dat1 = c(sig_diff1, sig_diff2, sig_diff3)
}else if ((length(language)==1 && length(output)>=9) | (length(language)==2 && length(output)>=8)){
dat1 = c(sig_diff1, sig_diff2)
}else{
dat1 = sig_diff1
}
assign(paste0("int_sig_diff_", var1, "_", var2), dat1, envir = .GlobalEnv)
rm(dat1, output)
# run plot_smooth to grab data for polar plots
if (condition=="tokenPooled.ord" && length(language)==1){
cond_p1 = capture.output(cat(paste0("list(", condition, "=var1)")))
cond_p2 = capture.output(cat(paste0("list(", condition, "=var2)")))
}else if (condition=="langNoteInt.ord" && length(language)==2){
cond_p1 = capture.output(cat(paste0("list(", condition, "='", language[1], ".", var1, ".", var2,
"', ", condition, "='", language[2], ".", var1, ".", var2, "')")))
cond_p2 = capture.output(cat(paste0("list(", condition, "='", language[2], ".", var1, ".", var2,
"', ", condition, "='", language[1], ".", var1, ".", var2, "')")))
}else if (condition=="native_lg.ord"){
cond_p1 = capture.output(cat(paste0("list(", condition, "=var1)")))
cond_p2 = capture.output(cat(paste0("list(", condition, "=var2)")))
}
# hardcoded condition
p1 = plot_smooth(x=get(model), view=values, cond=eval(parse(text=cond_p1)))
p2 = plot_smooth(x=get(model), view=values, cond=eval(parse(text=cond_p2)))
# set Rho max to the max of the fit + half the difference between max of the fit and the upper limit
maximum=max(p1$fv$fit, p2$fv$fit)+((max(p1$fv$ul, p2$fv$ul)-max(p1$fv$fit, p2$fv$fit))/2)
# plot in polar coordinates
if (exists("sig_diff4")){
p=plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff1[1]*180/pi, sig_diff1[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff2[1]*180/pi, sig_diff2[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff3[1]*180/pi, sig_diff3[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff4[1]*180/pi, sig_diff4[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)),
title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
Sys.sleep(0)
print(p)
}else if (exists("sig_diff3")){
p=plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff1[1]*180/pi, sig_diff1[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff2[1]*180/pi, sig_diff2[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff3[1]*180/pi, sig_diff3[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)),
title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
Sys.sleep(0)
print(p)
}else if (exists("sig_diff2")){
p=plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff1[1]*180/pi, sig_diff1[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff2[1]*180/pi, sig_diff2[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)),
title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
Sys.sleep(0)
print(p)
}else{
p=plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$fit, line=list(color="blue", width=2.5), name=names_smooths[1]) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ul, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p1$fv$theta_uncut_z*180/pi, r=p1$fv$ll, line=list(color="blue", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$fit, line=list(color="red", width=2.5), name=names_smooths[2]) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ul, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=p2$fv$theta_uncut_z*180/pi, r=p2$fv$ll, line=list(color="red", dash="dot", width=0.5), showlegend=FALSE) %>%
add_trace(theta=seq(sig_diff1[1]*180/pi, sig_diff1[2]*180/pi, length.out=20),
r=c(0, rep(maximum, 18), 0), line=list(color="black", width=0.5), fill="toself", fillcolor=rgb(0,0,0,max=255,alpha=25), showlegend=FALSE) %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)),
angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)),
title=plot_title, legend=list(orientation="h", xanchor="center", x=0.5))
Sys.sleep(0)
print(p)
}
}
}
df <- read.csv("all_data_NZE_Tongan_w_context_checked_7_March_2019_not_cut_checked_27_Apr.csv", sep=',', stringsAsFactors = F)
# remove empty column
df$X = NULL
df$tokenPooled <- factor(df$tokenPooled)
df$subject <- factor(df$subject)
df$sex <- factor(df$sex)
df$native_lg <- factor(df$native_lg)
# updated
df$playing_proficiency[df$playing_proficiency == "intermediate"] <- "amateur"
df$playing_proficiency <- factor(df$playing_proficiency, levels = c("amateur","semi-professional","professional"))
df$block <- factor(df$block)
df$theta_uncut_z <- as.numeric(df$theta_uncut_z)
# updated
df$note_intensity <- factor(df$note_intensity, levels = c("piano","mezzopiano","mezzoforte","forte"))
df$tokenPooled <- factor(df$tokenPooled)
dfSummary <- group_by(df, subject, playing_proficiency, activity, native_lg, tokenPooled, theta_uncut_z_group = cut(theta_uncut_z,breaks=100)) %>% summarise(rhoVar = sd(rho_uncut_z,na.rm=TRUE))
dfSummary$theta_uncut_z = as.character(dfSummary$theta_uncut_z_group)
dfSummary$theta_uncut_z = gsub("\\[|\\]|\\(|\\)", "",dfSummary$theta_uncut_z)
dfSummary$theta_uncut_z = strsplit(dfSummary$theta_uncut_z,",")
dfSummary$theta_uncut_z2 = 1
for(i in c(1:nrow(dfSummary)))
{
dfSummary$theta_uncut_z2[i] = mean(as.numeric(unlist(dfSummary$theta_uncut_z[i])))
}
dfSummary$theta_uncut_z=dfSummary$theta_uncut_z2
dfSummary$playing_proficiency.ord <- as.ordered(dfSummary$playing_proficiency)
contrasts(dfSummary$playing_proficiency.ord) <- "contr.treatment"
dfSummary$activity.ord <- as.ordered(dfSummary$activity)
contrasts(dfSummary$activity.ord) <- "contr.treatment"
dfSummary$langNoteInt <- interaction(dfSummary$native_lg,
dfSummary$tokenPooled)
dfSummary$langNoteInt.ord <- as.ordered(dfSummary$langNoteInt)
contrasts(dfSummary$langNoteInt.ord) <- "contr.treatment"
dfSummary$native_lg.ord <- as.ordered(dfSummary$native_lg)
contrasts(dfSummary$native_lg.ord) <- "contr.treatment"
dfSummary$tokenPooled.ord <- as.ordered(dfSummary$tokenPooled)
contrasts(dfSummary$tokenPooled.ord) <- "contr.treatment"
dfSummary = na.omit(dfSummary)
for(i in unique(dfSummary$subject))
{
dfSubject = subset(dfSummary,dfSummary$subject == i)
for(j in unique(dfSubject$tokenPooled))
{
dfSummary$start[dfSummary$subject == i & dfSummary$tokenPooled == j] <-
dfSummary$theta_uncut_z[dfSummary$subject == i & dfSummary$tokenPooled == j] ==
min(dfSummary$theta_uncut_z[dfSummary$subject == i & dfSummary$tokenPooled == j])
}
}
if (run_models == TRUE){
mdl.sys.time1 <- system.time(VAR.gam.noAR.Mod1 <- bam(rhoVar ~ langNoteInt.ord + s(theta_uncut_z, bs="cr", k=10) + s(theta_uncut_z, bs="cr", k=10, by=langNoteInt.ord) +s(theta_uncut_z, subject, bs="fs", k=10, m=1, by=langNoteInt.ord), data=dfSummary, discrete=TRUE, nthreads=ncores, method = "fREML"))
mdl.sys.time1
# save model & model summary so they can be reloaded later
saveRDS(VAR.gam.noAR.Mod1, paste0(output_dir,"/VAR.gam.noAR.Mod1.rds"))
capture.output(summary(VAR.gam.noAR.Mod1),
file = paste0(output_dir,"/summary_VAR.gam.noAR.Mod1.txt"))
}else{
# reload model from output_dir
VAR.gam.noAR.Mod1 = readRDS(paste0(output_dir,"/VAR.gam.noAR.Mod1.rds"))
}
if (run_models == TRUE){
# set start_value_rho from VAR.gam.noAR.Mod1 as rho_est for the next model
rho_est <- start_value_rho(VAR.gam.noAR.Mod1)
mdl.sys.time2 <- system.time(VAR.gam.AR.Mod2 <- bam(rhoVar ~ langNoteInt.ord + s(theta_uncut_z, bs="cr", k=10) + s(theta_uncut_z, bs="cr", k=10, by=langNoteInt.ord) + s(theta_uncut_z, subject, bs="fs", k=10, m=1, by=langNoteInt.ord), AR.start=dfSummary$start, rho=rho_est, data=dfSummary, discrete=TRUE, nthreads=ncores, method = "fREML"))
mdl.sys.time2
# save model & model summary so they can be reloaded later
saveRDS(VAR.gam.AR.Mod2, paste0(output_dir,"/VAR.gam.AR.Mod2.rds"))
capture.output(summary(VAR.gam.AR.Mod2),
file = paste0(output_dir,"/summary_VAR.gam.AR.Mod2.txt"))
}else{
# reload model from output_dir
VAR.gam.AR.Mod2 = readRDS(paste0(output_dir,"/VAR.gam.AR.Mod2.rds"))
}
acf_resid(VAR.gam.AR.Mod2, main = "Average ACF AR", cex.lab=1.5, cex.axis=1.5)
acf_resid(VAR.gam.AR.Mod2,split_pred=list(dfSummary$tokenPooled),main = "Average ACF AR by tokenPooled",cex.lab=1.5,cex.axis=1.5)
acf_resid(VAR.gam.AR.Mod2,split_pred=list(dfSummary$native_lg),main = "Average ACF AR by native language",cex.lab=1.5,cex.axis=1.5)
summary(VAR.gam.AR.Mod2)
Family: gaussian
Link function: identity
Formula:
rhoVar ~ langNoteInt.ord + s(theta_uncut_z, bs = "cr",
k = 10) + s(theta_uncut_z, bs = "cr", k = 10, by = langNoteInt.ord) +
s(theta_uncut_z, subject, bs = "fs", k = 10, m = 1,
by = langNoteInt.ord)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.00231 0.10777 74.256 < 2e-16 ***
langNoteInt.ordNZE.É\u0090 -0.97363 1.61302 -0.604 0.546111
langNoteInt.ordNZE.É\u0090Ë\u0090 -2.64035 1.85143 -1.426 0.153854
langNoteInt.ordNZE.É’ -2.46702 2.76454 -0.892 0.372202
langNoteInt.ordNZE.Bb2 -1.98406 1.29226 -1.535 0.124720
langNoteInt.ordTongan.Bb2 -1.41836 1.41220 -1.004 0.315220
langNoteInt.ordNZE.Bb3 -2.98062 1.73953 -1.713 0.086644 .
langNoteInt.ordTongan.Bb3 -3.58208 2.20014 -1.628 0.103519
langNoteInt.ordNZE.D4 -4.14783 1.66244 -2.495 0.012604 *
langNoteInt.ordTongan.D4 -5.16170 1.97641 -2.612 0.009018 **
langNoteInt.ordNZE.e 0.50229 1.63619 0.307 0.758855
langNoteInt.ordTongan.e -2.15181 0.92909 -2.316 0.020568 *
langNoteInt.ordNZE.É™ -0.07156 1.26678 -0.056 0.954950
langNoteInt.ordNZE.É™# -1.81077 0.40448 -4.477 7.63e-06 ***
langNoteInt.ordNZE.É› -1.25424 1.83884 -0.682 0.495196
langNoteInt.ordNZE.ɘ 0.04167 2.05228 0.020 0.983800
langNoteInt.ordNZE.F3 -3.09511 1.73763 -1.781 0.074894 .
langNoteInt.ordTongan.F3 -2.53651 1.38509 -1.831 0.067074 .
langNoteInt.ordNZE.F4 -5.29805 1.44440 -3.668 0.000245 ***
langNoteInt.ordTongan.F4 -8.89295 2.81838 -3.155 0.001606 **
langNoteInt.ordTongan.i -1.95829 1.33837 -1.463 0.143434
langNoteInt.ordNZE.iË\u0090 -29.28668 16.83868 -1.739 0.082009 .
langNoteInt.ordTongan.o -0.92331 5.02574 -0.184 0.854239
langNoteInt.ordNZE.oË\u0090 -2.13673 1.64411 -1.300 0.193747
langNoteInt.ordNZE.ɵË\u0090 -1.62651 2.37389 -0.685 0.493247
langNoteInt.ordTongan.u -6.83503 10.71965 -0.638 0.523732
langNoteInt.ordNZE.ʉË\u0090 -3.21586 3.80494 -0.845 0.398022
langNoteInt.ordNZE.ÊŠ -1.02949 3.74579 -0.275 0.783444
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df
s(theta_uncut_z) 8.872 8.936
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090 5.546 5.946
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090Ë\u0090 5.237 5.619
s(theta_uncut_z):langNoteInt.ordNZE.É’ 4.973 5.297
s(theta_uncut_z):langNoteInt.ordNZE.Bb2 4.453 4.756
s(theta_uncut_z):langNoteInt.ordTongan.Bb2 5.678 6.024
s(theta_uncut_z):langNoteInt.ordNZE.Bb3 1.003 1.004
s(theta_uncut_z):langNoteInt.ordTongan.Bb3 7.791 8.078
s(theta_uncut_z):langNoteInt.ordNZE.D4 2.924 3.092
s(theta_uncut_z):langNoteInt.ordTongan.D4 3.844 4.110
s(theta_uncut_z):langNoteInt.ordNZE.e 6.911 7.197
s(theta_uncut_z):langNoteInt.ordTongan.e 6.022 6.443
s(theta_uncut_z):langNoteInt.ordNZE.É™ 8.038 8.378
s(theta_uncut_z):langNoteInt.ordNZE.É™# 8.412 8.698
s(theta_uncut_z):langNoteInt.ordNZE.É› 6.415 6.796
s(theta_uncut_z):langNoteInt.ordNZE.ɘ 7.059 7.409
s(theta_uncut_z):langNoteInt.ordNZE.F3 1.001 1.001
s(theta_uncut_z):langNoteInt.ordTongan.F3 6.412 6.792
s(theta_uncut_z):langNoteInt.ordNZE.F4 4.088 4.346
s(theta_uncut_z):langNoteInt.ordTongan.F4 4.886 5.195
s(theta_uncut_z):langNoteInt.ordTongan.i 7.799 8.150
s(theta_uncut_z):langNoteInt.ordNZE.iË\u0090 6.118 6.552
s(theta_uncut_z):langNoteInt.ordTongan.o 6.065 6.458
s(theta_uncut_z):langNoteInt.ordNZE.oË\u0090 5.824 6.194
s(theta_uncut_z):langNoteInt.ordNZE.ɵË\u0090 1.001 1.001
s(theta_uncut_z):langNoteInt.ordTongan.u 6.249 6.720
s(theta_uncut_z):langNoteInt.ordNZE.ʉË\u0090 1.001 1.001
s(theta_uncut_z):langNoteInt.ordNZE.ÊŠ 2.316 2.436
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090 70.857 100.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090Ë\u0090 70.884 100.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É’ 79.953 100.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb2 66.321 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb2 66.393 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb3 71.568 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb3 69.313 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.D4 67.242 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.D4 68.879 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.e 66.160 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.e 61.414 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™ 59.777 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™# 48.910 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.É› 69.393 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɘ 71.632 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.F3 75.905 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.F3 66.936 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.F4 61.214 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.F4 66.915 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.i 62.482 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.iË\u0090 79.667 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.o 82.606 100.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.oË\u0090 68.673 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɵË\u0090 77.097 99.000
s(theta_uncut_z,subject):langNoteInt.ordTongan.u 84.433 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.ʉË\u0090 81.023 99.000
s(theta_uncut_z,subject):langNoteInt.ordNZE.ÊŠ 77.375 99.000
F p-value
s(theta_uncut_z) 70.489 < 2e-16
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090 1.445 0.169930
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090Ë\u0090 2.493 0.025088
s(theta_uncut_z):langNoteInt.ordNZE.É’ 1.250 0.205920
s(theta_uncut_z):langNoteInt.ordNZE.Bb2 1.584 0.301922
s(theta_uncut_z):langNoteInt.ordTongan.Bb2 2.210 0.041133
s(theta_uncut_z):langNoteInt.ordNZE.Bb3 0.093 0.762294
s(theta_uncut_z):langNoteInt.ordTongan.Bb3 3.113 0.000804
s(theta_uncut_z):langNoteInt.ordNZE.D4 1.113 0.345495
s(theta_uncut_z):langNoteInt.ordTongan.D4 1.625 0.157326
s(theta_uncut_z):langNoteInt.ordNZE.e 3.014 0.004660
s(theta_uncut_z):langNoteInt.ordTongan.e 1.728 0.073853
s(theta_uncut_z):langNoteInt.ordNZE.É™ 8.910 1.75e-11
s(theta_uncut_z):langNoteInt.ordNZE.É™# 18.668 < 2e-16
s(theta_uncut_z):langNoteInt.ordNZE.É› 2.142 0.067326
s(theta_uncut_z):langNoteInt.ordNZE.ɘ 3.775 0.000484
s(theta_uncut_z):langNoteInt.ordNZE.F3 0.057 0.812421
s(theta_uncut_z):langNoteInt.ordTongan.F3 2.265 0.018740
s(theta_uncut_z):langNoteInt.ordNZE.F4 1.708 0.170357
s(theta_uncut_z):langNoteInt.ordTongan.F4 2.669 0.033058
s(theta_uncut_z):langNoteInt.ordTongan.i 3.002 0.001150
s(theta_uncut_z):langNoteInt.ordNZE.iË\u0090 2.756 0.006887
s(theta_uncut_z):langNoteInt.ordTongan.o 2.203 0.048544
s(theta_uncut_z):langNoteInt.ordNZE.oË\u0090 2.666 0.012507
s(theta_uncut_z):langNoteInt.ordNZE.ɵË\u0090 0.026 0.873668
s(theta_uncut_z):langNoteInt.ordTongan.u 2.042 0.055998
s(theta_uncut_z):langNoteInt.ordNZE.ʉË\u0090 0.413 0.520613
s(theta_uncut_z):langNoteInt.ordNZE.ÊŠ 0.180 0.783067
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090 6.865 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090Ë\u0090 7.818 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É’ 18.664 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb2 10.303 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb2 18.281 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb3 14.466 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb3 22.632 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.D4 13.933 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.D4 25.480 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.e 12.291 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.e 8.573 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™ 14.365 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™# 3.117 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.É› 8.070 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɘ 13.504 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.F3 15.014 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.F3 20.935 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.F4 11.949 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.F4 26.468 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.i 14.004 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.iË\u0090 45.669 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.o 22.896 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.oË\u0090 13.644 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɵË\u0090 12.675 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordTongan.u 43.466 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.ʉË\u0090 22.778 < 2e-16
s(theta_uncut_z,subject):langNoteInt.ordNZE.ÊŠ 28.619 < 2e-16
s(theta_uncut_z) ***
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090
s(theta_uncut_z):langNoteInt.ordNZE.É\u0090Ë\u0090 *
s(theta_uncut_z):langNoteInt.ordNZE.É’
s(theta_uncut_z):langNoteInt.ordNZE.Bb2
s(theta_uncut_z):langNoteInt.ordTongan.Bb2 *
s(theta_uncut_z):langNoteInt.ordNZE.Bb3
s(theta_uncut_z):langNoteInt.ordTongan.Bb3 ***
s(theta_uncut_z):langNoteInt.ordNZE.D4
s(theta_uncut_z):langNoteInt.ordTongan.D4
s(theta_uncut_z):langNoteInt.ordNZE.e **
s(theta_uncut_z):langNoteInt.ordTongan.e .
s(theta_uncut_z):langNoteInt.ordNZE.É™ ***
s(theta_uncut_z):langNoteInt.ordNZE.É™# ***
s(theta_uncut_z):langNoteInt.ordNZE.É› .
s(theta_uncut_z):langNoteInt.ordNZE.ɘ ***
s(theta_uncut_z):langNoteInt.ordNZE.F3
s(theta_uncut_z):langNoteInt.ordTongan.F3 *
s(theta_uncut_z):langNoteInt.ordNZE.F4
s(theta_uncut_z):langNoteInt.ordTongan.F4 *
s(theta_uncut_z):langNoteInt.ordTongan.i **
s(theta_uncut_z):langNoteInt.ordNZE.iË\u0090 **
s(theta_uncut_z):langNoteInt.ordTongan.o *
s(theta_uncut_z):langNoteInt.ordNZE.oË\u0090 *
s(theta_uncut_z):langNoteInt.ordNZE.ɵË\u0090
s(theta_uncut_z):langNoteInt.ordTongan.u .
s(theta_uncut_z):langNoteInt.ordNZE.ʉË\u0090
s(theta_uncut_z):langNoteInt.ordNZE.ÊŠ
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É\u0090Ë\u0090 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É’ ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb2 ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb2 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.Bb3 ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.Bb3 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.D4 ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.D4 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.e ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.e ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™ ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É™# ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.É› ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɘ ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.F3 ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.F3 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.F4 ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.F4 ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.i ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.iË\u0090 ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.o ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.oË\u0090 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.ɵË\u0090 ***
s(theta_uncut_z,subject):langNoteInt.ordTongan.u ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.ʉË\u0090 ***
s(theta_uncut_z,subject):langNoteInt.ordNZE.ÊŠ ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.868 Deviance explained = 88.2%
fREML = 36912 Scale est. = 2.0136 n = 19150
And finally overall variance differences by language
# make copy of dfNotes
dat1 = dfSummary
dat1$predicted_values = predict(VAR.gam.AR.Mod2)
# plot in polar coordinates using plotly
dat1_NZE = dat1[dat1$native_lg == "NZE",]
dat1_Tongan = dat1[dat1$native_lg == "Tongan",]
# estimate smooths using R's generic predict.smooth.spline function
smooth_NZE=list(theta=seq(min(dat1_NZE$theta_uncut_z)*180/pi, max(dat1_NZE$theta_uncut_z)*180/pi, length=100), r=predict(smooth.spline(dat1_NZE$theta_uncut_z, dat1_NZE$predicted_values),
seq(min(dat1_NZE$theta_uncut_z), max(dat1_NZE$theta_uncut_z), length=100))$y, line=list(color="blue", dash="dash"))
smooth_Tongan=list(theta=seq(min(dat1_Tongan$theta_uncut_z)*180/pi, max(dat1_Tongan$theta_uncut_z)*180/pi, length=100), r=predict(smooth.spline(dat1_Tongan$theta_uncut_z, dat1_Tongan$predicted_values),
seq(min(dat1_Tongan$theta_uncut_z), max(dat1_Tongan$theta_uncut_z), length=100))$y, line=list(color="red", dash=""))
# set Rho max to the max of the predicted_values + 5
max_predictions = max(max(smooth_NZE$r), max(smooth_Tongan$r))
maximum=max_predictions+5
rm(max_predictions)
p = plot_ly(type='scatterpolar', mode='lines') %>%
add_trace(theta=smooth_NZE$theta, r=smooth_NZE$r, line=list(color=smooth_NZE$line$color[[1]], width=2.5, dash=smooth_NZE$line$dash[[1]]), name="overall average of NZE notes") %>%
add_trace(theta=smooth_Tongan$theta, r=smooth_Tongan$r, line=list(color=smooth_Tongan$line$color[[1]], width=2.5, dash=smooth_Tongan$line$dash[[1]]), name="overall average of Tongan notes") %>%
layout(polar=list(sector=c(20,160), radialaxis=list(angle=90, range=c(0,maximum)), angularaxis=list(thetaunit='radians', direction="clockwise", rotation=0)), title="Overall average smooths for NZE vs Tongan note productions", legend=list(orientation="h", xanchor="center", x=0.5))
p
rm(dat1, dat1_NZE, dat1_Tongan, maximum, smooth_NZE, smooth_Tongan)